{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Solving Non-linear Equations\n",
    "\n",
    "The [Roots](https://github.com/JuliaMath/Roots.jl) package provides methods for solving a non-linear equation (one variable, one function). \n",
    "\n",
    "For a system of non-linear equations, use [NLsolve](https://github.com/JuliaNLSolvers/NLsolve.jl)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "printyellow (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Dates, Roots, NLsolve\n",
    "\n",
    "include(\"printmat.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "#pyplot(size=(600,400))\n",
    "gr(size=(480,320))\n",
    "default(fmt = :svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining and Plotting the Function\n",
    "\n",
    "The next few cells define a fairly simple function and then plots it. \n",
    "\n",
    "If possible, plot your function. Maybe you see something strange, perhaps there are several roots? It also helps you set the initial guesses (or brackets) for root solving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fn1 (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fn1(x,c) = 2*(x - 1.1)^2 - c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"480\" height=\"320\" viewBox=\"0 0 1920 1280\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip6200\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"1920\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6200)\" d=\"\n",
       "M0 1280 L1920 1280 L1920 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6201\">\n",
       "    <rect x=\"384\" y=\"0\" width=\"1345\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6200)\" d=\"\n",
       "M175.611 1105.62 L1872.76 1105.62 L1872.76 121.675 L175.611 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6202\">\n",
       "    <rect x=\"175\" y=\"121\" width=\"1698\" height=\"985\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  223.643,1105.62 223.643,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  623.913,1105.62 623.913,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1024.18,1105.62 1024.18,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1424.45,1105.62 1424.45,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1824.72,1105.62 1824.72,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  175.611,1025.15 1872.76,1025.15 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  175.611,814.662 1872.76,814.662 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  175.611,604.175 1872.76,604.175 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  175.611,393.688 1872.76,393.688 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  175.611,183.201 1872.76,183.201 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,1105.62 1872.76,1105.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,1105.62 175.611,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  223.643,1105.62 223.643,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  623.913,1105.62 623.913,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1024.18,1105.62 1024.18,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1424.45,1105.62 1424.45,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1824.72,1105.62 1824.72,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,1025.15 201.068,1025.15 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,814.662 201.068,814.662 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,604.175 201.068,604.175 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,393.688 201.068,393.688 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  175.611,183.201 201.068,183.201 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 223.643, 1157.22)\" x=\"223.643\" y=\"1157.22\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 623.913, 1157.22)\" x=\"623.913\" y=\"1157.22\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1024.18, 1157.22)\" x=\"1024.18\" y=\"1157.22\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1424.45, 1157.22)\" x=\"1424.45\" y=\"1157.22\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1824.72, 1157.22)\" x=\"1824.72\" y=\"1157.22\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 156.411, 1042.65)\" x=\"156.411\" y=\"1042.65\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 156.411, 832.162)\" x=\"156.411\" y=\"832.162\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 156.411, 621.675)\" x=\"156.411\" y=\"621.675\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 156.411, 411.188)\" x=\"156.411\" y=\"411.188\">6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 156.411, 200.701)\" x=\"156.411\" y=\"200.701\">8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1024.18, 73.2)\" x=\"1024.18\" y=\"73.2\">the fn1(x,0.5) function</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1024.18, 1234.68)\" x=\"1024.18\" y=\"1234.68\">x</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 98.8861, 613.647)\" x=\"98.8861\" y=\"613.647\">y</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip6202)\" style=\"stroke:#ff0000; stroke-width:8; stroke-opacity:1; fill:none\" points=\"\n",
       "  223.643,149.523 263.67,235.822 303.697,317.912 343.724,395.793 383.751,469.463 423.778,538.924 463.805,604.175 503.832,665.216 543.859,722.048 583.886,774.67 \n",
       "  623.913,823.082 663.94,867.284 703.967,907.277 743.994,943.059 784.021,974.632 824.048,1002 864.075,1025.15 904.102,1044.09 944.129,1058.83 984.156,1069.35 \n",
       "  1024.18,1075.67 1064.21,1077.77 1104.24,1075.67 1144.26,1069.35 1184.29,1058.83 1224.32,1044.09 1264.35,1025.15 1304.37,1002 1344.4,974.632 1384.43,943.059 \n",
       "  1424.45,907.277 1464.48,867.284 1504.51,823.082 1544.53,774.67 1584.56,722.048 1624.59,665.216 1664.62,604.175 1704.64,538.924 1744.67,469.463 1784.7,395.793 \n",
       "  1824.72,317.912 \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = -1:0.1:3\n",
    "\n",
    "plot( x,fn1.(x,0.5),\n",
    "      linecolor = :red,\n",
    "      linewidth = 2,\n",
    "      legend = nothing,\n",
    "      title = \"the fn1(x,0.5) function\",\n",
    "      xlabel = \"x\",\n",
    "      ylabel = \"y\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There seems to be two roots: around 0.6 and 1.6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a Non-Linear Equation\n",
    "\n",
    "The Roots package wants a function with only one input. An easy way to turn ```fn1(a,0.5)``` into that form is by defining an anonymous function:\n",
    "```\n",
    "x->fn1(x,0.5)\n",
    "```\n",
    "\n",
    "Then, running \n",
    "```\n",
    "find_zero(x->fn1(x,0.5),(x₀,x₁))\n",
    "```\n",
    "searches for a root in the `[x₀,x₁]` interval. Alternatively, you can also do \n",
    "```\n",
    "find_zero(x->fn1(x,0.5),x₂)\n",
    "``` \n",
    "where `x₂` is a single starting guess.\n",
    "\n",
    "Instead, running\n",
    "```\n",
    "find_zeros(x->fn1(x,0.5),x₀,x₁)\n",
    "```\n",
    "searches for all roots between x₀ and x₁. (Notice the *s* in `find_zeros`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "at which x is fn1(x,0.5) = 0?      0.600\n",
      "at which x is fn1(x,0.5) = 0?      1.600\n",
      "\n",
      "\u001b[34m\u001b[1myes, there are several roots. Just look at it (in the plot)\u001b[22m\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "xRoot1 = find_zero(x->fn1(x,0.5),(-1,1))            #searches for roots in [-1,1]\n",
    "printlnPs(\"at which x is fn1(x,0.5) = 0? \",xRoot1)\n",
    "\n",
    "xRoot2 = find_zero(x->fn1(x,0.5),2)              #searches for roots around 2\n",
    "printlnPs(\"at which x is fn1(x,0.5) = 0? \",xRoot2)\n",
    "\n",
    "printblue(\"\\nyes, there are several roots. Just look at it (in the plot)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "at which x is fn1(x,0.5) = 0?      0.600     1.600\n"
     ]
    }
   ],
   "source": [
    "xRootsAll = find_zeros(x->fn1(x,0.5),-1,3)            #find_zeros (notice the \"s\")\n",
    "\n",
    "printlnPs(\"at which x is fn1(x,0.5) = 0? \",xRootsAll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a System of Non-Linear Equations\n",
    "\n",
    "The NLsolve package has many options. The cells below illustrate a very simple case (2 non-linear equations with 2 unknowns, no information about the derivatives).\n",
    "\n",
    "The two equations are\n",
    "\n",
    "$ \n",
    "y-x^2-1=0\n",
    "$\n",
    "\n",
    "$\n",
    "y-x-1=0\n",
    "$\n",
    "\n",
    "and the roots are at $(x,y)=(0,1)$ and also at $(1,2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fn2 (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fn2(p)             #p is a vector with 2 elements, the output too\n",
    "    (x,y) = (p[1],p[2])\n",
    "    z = [y-x^2-1;y-x-1]     #equal to [0,0]\n",
    "    return z\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There is a solution at             -0.000     1.000\n",
      "There is a another solution at      1.000     2.000\n"
     ]
    }
   ],
   "source": [
    "Sol2a = nlsolve(fn2,[0.0,0.5])\n",
    "printlnPs(\"There is a solution at         \",Sol2a.zero)\n",
    "\n",
    "Sol2b = nlsolve(fn2,[1.0,0.0])  #try again, using another starting guess\n",
    "printlnPs(\"There is a another solution at \",Sol2b.zero)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.3.0",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
